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Weighted scale-free networks with topology-dependent interactions are studied. It is shown that 
the possible universality classes of critical behaviour, which are known to depend on topology, can 
also be explored by tuning the form of the interactions at fixed topology. For a model of opinion 
formation, simple mean field and scaling arguments show that a mapping 7' = (7 — /n)/(l — fi) 
. describes how a shift of the standard exponent 7 of the degree distribution can absorb the effect 

1) ' of degree-dependent pair interactions Jij <x (fc;fcj) _M , where hi stands for the degree of vertex i. 

This prediction is verified by extensive numerical investigations using the cavity method and Monte 
Carlo simulations. The critical temperature of the model is obtained through the Bethe-Peierls 
approximation and with the replica technique. The mapping can be extended to nonequilibrium 
models such as those describing the spreading of a disease on a network. 
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I. INTRODUCTION 



In recent years 0,|3it has become clear that many natural and technological networks are quite different from simple 
random networks [Sand share several unexpected properties such as a scale-free degree distribution 0, 0, small 
world connectivity |4| , soft modularity^?] and so on. Much work has gone into a precise characterisation of these 
topological properties for a variety of networks as diverse as the internet |fj , metabolic networks in cells or networks 
of chemical reactions in planetary atmospheres Q . Nontrivial topology is now established as an essential ingredient of 
\Q , complex systems |9j. This observation naturally raises the question how these structures are formed and grow. A well 
known mechanism is that of preferential attachment which leads to the Barabasi- Albert network with a power-law 
degree distribution 

Another important question is how the network topology affects physical properties such as collective behaviour, 
transport quantities, the spreading of a disturbance, ... Particular interest has been devoted to the behaviour of 
the Ising model on a network. Besides being the standard model of (equilibrium) statistical mechanics, the Ising 

S model is also expected to give a simple description of sociological phenomena such as opinion formation . The 
first studies of this model indicated pjj that on a Barabasi- Albert network, the Ising model is always ordered. It 
was, however, soon realised that a finite transition temperature can be obtained if one considers the Ising model on 
scale- free networks with a degree distribution that is less long ranged than that of the Barabasi- Albert one [iH . Il4| . 
Even more interesting is the observation of non-trivial critical behaviour for these cases [l3l Il4l ITsL Il6| . 

Another mechanism that leads to a finite transition temperature was discovered by one of us in the so called special 
attention network. This is a model of a weighted network where the interactions between connected spins are made 
topology-dependent ^7|- A further investigation of that and related models then led to the discovery that topology 
and interaction can be traded', in the sense that the effect of a change in interaction law can be transformed away by 
an appropriate change of the degree distribution law. In the present paper we present a detailed investigation of the 
critical behaviour of the Ising model with degree-dependent interactions. We also discuss the extension of our results 
to nonequilibrium situations. A summary of our results was published earlier |18| . 

This paper is organised as follows. In the next section, we define the model and present the results of a simple mean 
field theory. In section 3, we study the model with the cavity approach. In section 4, we describe the application of the 
Bethe-Peierls method and the replica technique to our model. In section 5, we present the results of extensive Monte 
Carlo simulations. In section 6, we discuss the extension of our main result to a nonequilibrium process. Finally, we 
present our conclusions in section 7. 



2 



II. THE MODEL 



Our model can be defined on a general network (or graph). When in a graph a vertex (or node) i is connected with 
ki other nodes we say it has degree fcj. We will mostly have in mind a scale-free network for which the degree is a 
random variable whose distribution P(k) is a power-law: P(k) oc fc~ 7 . We take 7 > 2 so that the average degree 
Q = J kP(k)dk is finite. For future reference we also introduce the notation Qi for the second moment of P(k). The 
Barabasi- Albert (BA) network has 7 = 3 10] . The number of nodes in the network is denoted by N. 

We next define an Ising model on this network by associating to each node i a variable Sj = ±1 and to each pair 
of linked nodes an energy —JijSiSj. In this paper, we will choose the couplings to be given by 

Jn = JQ^/ihkjY ■ (l) 

with J > 0. For /1 = 0, the behaviour of this model has been investigated by various authors and the critical behaviour 
was found to depend on 7 Q E [j^ 

For 7 > 5 the results of standard mean field theory were found to apply: the 
critical temperature T c is finite, the order parameter vanishes with an exponent 1/2 and the specific heat has a finite 
jump at T c . When 3 < 7 < 5, the transition temperature remains finite, the specific heat goes to zero continuously 
as a function of temperature with an exponent that depends on 7. The same is true for the order parameter. In the 
borderline case 7 = 5, logarithmic corrections appear. For all cases with 7 > 3, the zero-field susceptibility diverges 
as \T — T c | _1 . Finally, when 2 < 7 < 3, the system is ordered at all temperatures. 

When [i = 1/2, our model corresponds to the special attention network (SAN) introduced earlier by one of us [l7| . 
A simple mean field approach showed that on a BA network, the SAN has a finite transition temperature 0]. We 
start by extending that simple approach to the case of general /i and 7. 



A. A simple mean-field approximation 



For a given network realization, the local magnetisation m, = (sj) at node i obeys the exact equation 

( 



\3 = 1 



Here (.) is the thermal average, ks is Boltzmann's constant and T is the temperature. Following Bianconi [ijj, we 
now apply a double mean field approximation to this equation by rewriting it as 



tanh [ y M m ] 



(3) 



where the sum now runs over all the nodes and [.] denotes the average over all the realisations of the network with a 
fixed set of degrees {ki}. For a BA network the probability pij that two nodes are connected was shown to be equal 
to kikj I (QN), 0] and we can expect this result to hold in some other cases ^(J. We therefore have 

[Jij] = -I.,!', J = (kikj) 1 ^ , (4) 

so that Eq. © can be rewritten as 

m, = tanh | — — — k) ^— k] ''to, I . (5) 

I k B T N ^ 3 ' 

The quantity 

OM-l N 

s =V^^' (6) 

J'=l 

is a convenient order parameter. From Eq. (J5J it follows that S obeys the selfconsistency equation 
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For N — > oo, the term on the right hand side is simply related to the average over the distribution P{k) so that we 
can rewrite Eq. as 



S = Q"- 1 [ fc^tanh ( 'I^-k^^s] P(k)dk 
im \ k BT J 



(8) 



Here m > 1 is the lowest degree that is possible in the network. Eq. JS} can be analysed in a standard way. For 
example, the critical temperature is determined by assuming S to be small. After linearisation we then obtain 

Tc = jq ^~ x r k 2 - 2 »p(k)dk . (9) 



For the power-law distribution this gives a finite T c provided 7 > 3 — 2fi. This is consistent with the earlier finding of 
a finite T c for the SAN on the BA network. One could then go on and determine, for example, the exponent /3 from 
a further analysis of Eq. ft is however more suitable to perform the transformation of variables 

k' = Q^k 1 -^ (10) 

This gives for the case that the degree distribution is power-law, and for fi < 1 

S = A (fc'^tanh (±— 1 )dk' (11) 
Jm' V k sT J 

where A is a constant and m! = Q^m}~^. Comparison with JSJ teaches us that (|ll|l is precisely the mean field 
equation for fi' — and for a degree distribution with an exponent that is modified to 

This relation can be expected to hold more generally. Indeed, it should be valid whenever within a mean-field approach 
the degree k only enters in physical properties through the quenched average interaction between any two nodes with 
fixed degrees, ki and kj. This average is, as shown above, JijPij with pij — kikj /(QN). For fi ^ 0, the ki can then be 
transformed using l|10|l . In order to retain the same physics, averages over the degree distribution must be invariant. 
This requires a distribution transformation 

P{k )=P'(k'(k))^- (13) 

from which Ijl2|l follows for a scale-free P(k). In section 6, we will in fact show that (|12|l also holds for a contact process 
with degree-dependent infection rates. The general derivation scheme outlined above, will also apply there. Thus, 
our simple mean field analysis indicates that by tuning /i G [2 — 7, 1] we can encounter the whole range of universality 
classes that was found in earlier work at (i = 0. As an example, one can study the whole set of universality classes by 
working on a BA-network and tuning fi. For the particular example of the SAN, l|12l) gives 7' = 5. 

In the rest of this paper we will apply several other approaches that go beyond simple mean field to our model. 
These will allow us to get a better estimate than (|5J) for the critical temperature, and to verify the equivalence between 
universality classes as expressed in Eq. |JT3J. 

We also observe that from Q and © it follows that 

"•■ = to ° "(^"^"E^ (14) 

where the last approximation holds close to T c . We therefore anticipate that the local magnetisation is proportional 
to k; ^ , a result that we will use further in this paper and which will be verified numerically. 



III. THE CAVITY METHOD 



We begin with a study of our model using the cavity approach. One reason for this is that concepts like the cavity 
field and the propagated field, which also appear in the replica study of our model, can be introduced in a physically 
transparent way within the context of this method. 
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The cavity method is closely related to the well known Bethe-Peierls (BP) approximation for a spin model on 
a Cay ley tree |20|. In fact, both methods are equivalent when there is replica symmetry (for the precise relation 
between various mean- field approximations, see reference 20). The cavity method was applied to diluted spin models 
in references 21-23. For Ising like models on a network, the method has the advantage that it can take into account 
site degree correlations, which certainly are present in networks grown via a preferential attachment rule pil |25| or 
in real world networks (from biology, sociology, information technology, ...). These correlations can be measured in 
terms of the clustering coefficient, the betweenness 1] or the abundance of cycles of a given length [2(|[27|]. This has 
to be contrasted with the replica approach of the next section, which, in its simplest form assumes independence of 
the site degrees. Extensions of the replica approach that take into account degree correlations are known [2q. but 
necessarily are more involved. Another attractive property of the cavity method is that thermodynamic quantities 
like the magnetisation can be calculated on a single network realisation. Together these properties allow the analysis 
of real world networks and to work with ensembles that, as in the BA case, are defined only through a growth process. 

To derive the cavity equations, it is assumed that the network has the structure of a tree. It is at first sight 
paradoxical that the tree approximation is adequate for determining the critical point of the network, because the 
Ising model on a tree (without loops) cannot display spontaneous symmetry breaking (SSB) at finite temperature 
and thus T c — 0. However, in the same way that the Bethe-Peierls approximation for an Ising model on a Cayley 
tree introduces an effective symmetry breaking field in the bulk so that SSB becomes possible, the cavity method 
introduces a random distribution of effective fields in the bulk (all fields being of the same sign) so that SSB becomes 
possible notwithstanding the absence of loops on the tree. Thus, the SSB due to the sparse loops in the actual network 
is replaced, in the cavity method, by the SSB due to the effective fields. The question remains whether the two SSB 
mechanisms lead to exactly the same value of T c . 

Consider then a particular site j in the network (see Fig. The assumption of tree structure implies that the 
sites connected to j form the roots of a set of independent subgraphs. 



The site j is connected to kj others. In the presence of an external field H the magnetisation at site j is given by 
(with f3 = l/k B T) 



Here Zij(s) is the partition sum of the whole subgraph starting from site i, and s is the value of the spin at j. This 
partition sum also includes the bond between the vertices i and j. Clearly, we can always write 




FIG. 1: Local structure of a network used in the derivation of the cavity equations (see text). 



m 



e m n*UZ l3 (+)-e-P H U*UZ^(-) 
e^ltiUZ^ + e-^ltfuZ^-) 



(15) 



Z i:j (s) «e^ s 

where we call uy the local cavity message. With this assumption, Eq. l|15f) becomes 



(16) 




(17) 



which also gives a clear physical interpretation to tiy. Indeed, Uij can be seen as the contribution to the total magnetic 
field acting on site j coming from the independent subtree having site i as root. In the cavity ansatz, those local 
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contributions are indeed independent. We need extra equations to determine the cavity fields selfconsistently. In 
order to determine these, we write Z^(s) in terms of the subtrees that are connected to it (see Fig. QJ. We denote 
the sites adjacent to i by the index I. Their number equals fej, but clearly, one of them is the starting site j. From 
the definition of Zij(s), we have 

^(») = ^'jj 1 Zl) ( + ) + e -Wvjj 1 ZK (_) (18) 
1=1 i=i 
so that using Eq. 115|l and Eq. (|16|l . we obtain after a little algebra 
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ij — — tanh (tanh (/3/i,y) tanh (/3Jy)) (19) 



where 



ki-l 

hio = H+j2u H (20) 
i=i 

is called the propagated field (sometimes also called cavity field). The equations (|19|) - <|20l) are known as belief 
propagation equations |20|. They were derived iteratively here, but it is possible to obtain the same equations 
through a variational procedure and assuming that configuration probabilities properly factorize on a tree. A short 
derivation from the variational approach will be presented in section 4.1. 

We solved the belief propagation equations iteratively for different values of /j, in our model. One can show in this 
case that belief propagation equations, when converging, do so to a unique fixed point (beside the fully paramagnetic 
solution which is always present in absence of external field). When possible, we could also be interested in averaging 
belief propagation equations over a proper random network ensemble. In that c ase, equations for propagated messages 
and fields become a set of integral equations for field probability distributions [22L l23| . This unfortunately cannot be 
easily done in the case of growing network ensembles, like the BA one. In the case of uncorrelated power-law degree 
distributed graphs, the self consistent integral equations can be written straightforwardly from (|19|) and (|2(Jfl . The 
resulting equations turn out to be exactly equations (|48|l and l|49|l that we will obtain through the replica approach. On 
the same random network ensemble, the replica symmetric calculation and the cavity method are therefore completely 
equivalent. Note, however, that the analytical results for T c presented in the replica calculation are valid for a graph 
ensemble which is not exactly the BA one. It is therefore plausible that we obtain a slight discrepancy between the 
replica and the cavity results obtained after averaging over real BA network realisations. 

Within the cavity approach, it is also possible to obtain the free energy F{(3) and from this other thermodynamic 
quantities such as the energy U(j3) and the specific heat C(/3). Here we only quote the result for U 

r/fm = _V7 ( t&n Hf3Jij) +tenHPh ij )t&nh(f3h ji ) \ 
[P> ^ lJ \l + tanh(/3Jy)tanh(/3/iy)tanh(^%)y 1 1 

We now turn to a discussion of our numerical results. These were obtained on BA networks. Ensemble averages 
in the case of BA networks can be performed numerically generating a large number of graphs with a given number 
of nodes with the usual preferential attachment rule, and subsequently averaging the results over all the graphs at 
a given N. Quantities such as T c and the critical exponents can then be found using finite size scaling. We studied 
1000 realisations of networks of N = 50, 100, 250, 500, 1000, 5000 and 10 4 nodes, 100 realisations with N = 5 x 10 4 
nodes and 10 realisations of 10 5 and 10 6 nodes. For each of these the cavity messages were determined iteratively. 
The numerical results showed little fluctuations when comparing different realisations, even for small network sizes. 

We first investigated the SAN and studied cases with Q = 4, 6, 8, 10 and 20. In table 1 (see section 4) we give our 
estimates for T c as a function of Q. These results will be compared with those coming from other methods. 

From the behaviour of the magnetisation M [M = -k J^i m i) below T c (Fig. we can determine the exponent j3. 
We find a value which is very close to the mean field value of 1/2. This value hardly depends on N or on Q. These 
data also allow us to obtain finite size estimates T C (N) for the critical temperature. We measured the magnetisation 
within a small temperature window around each value of T C (N) and this for small values of the external field H. 
From an extrapolation for large N , the value of the exponent 5 was found to be close to 3. More precisely we find 
1/5 = 0.333 ±0.010 for Q = 4. For bigger Q- values, we find the same value for 1/6 while the error decreases with Q. 
Using the scaling relation j s — [3(5 — 1), this leads to the value 1 for the susceptibility exponent 7 S . 

Specific heat computations (see Fig. |3J) show a jump around the transition. However, the discontinuity seems to 
become smaller for large N such that 

dC 

lim -== = -oo (22) 

T— >T~ U-£ 
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FIG. 2: Magnetisation versus temperature for Q = 4 (left) and Q = 10 (right) for AT = 10 2 , 10 3 , 10 4 , 10 5 . For Q = 4, results for 
N = 10 6 are also included. 



Both cases (vanishing jump with diverging derivative or discontinuity in C) can be consistent with a critical exponent 
a = 0. 

The conjecture embodied in H12[l predicts that the SAN on a BA network has 7' = 5. The critical exponents are 
then given by 7 S = 1, /? = 1/2 and a = [T3L flil HEj. In conclusion, we can say that the results from the cavity 
method for the SAN are in agreement with the relation (|12fl . 




FIG. 3: Specific heat versus temperature for Q = 4 (left) and Q = 10 (right) for N = 10 2 , 10 3 , 10 4 , 10 5 . For Q = 4, results for 
AT = 10 6 are also included. 

A peculiar behaviour of the specific heat can be seen in Fig. |3| For both Q~ vaiues shown there, one observes 
a maximum in the specific heat at a temperature below the critical one. This seems not to be a finite size effect. 
Moreover, as Fig. 0] shows, this kind of behaviour does not appear for the SAN on a network with a Poisson degree 
distribution. 

Another interesting feature is that close to T c we observe an inversion phenomenon in the local magnetisation of 
sites. Below the transition, there is a clear hierarchy in sites according to their spin magnetisation. Even though 
couplings are weighted, hubs are still the most magnetized sites and are believed to drive the transition. Plots of the 
local magnetization versus site degree for a few values of the temperature above and below T c are shown in Fig. [S] 
These data are consistent with the prediction (|14|) coming from the simple mean field theory. 

As can be seen in the scatter plots of Fig. [U low degree nodes connected to hubs have a lower magnetisation than 
equal degree nodes that are not connected to hubs. Above the transition, however, this situation is reversed. This is 
different from the situation for an Ising model with degree-independent couplings on a network, where at fixed degree, 
nodes connected to hubs are always most magnetised. 

Besides the SAN, we also investigated our model with the same techniques for // = 1/3 on a BA network. For that 
case, 111 'it predicts 7' = 4 which leads to the exponent values j s — (3 — —a = 1 and 6 = 2. In Fig. 0we show some 
representative results for the specific heat and the magnetisation as a function of temperature (both at Q = 4) . They 
show the appearance with increasing N of a regime where both quantities depend linearly on temperature, consistent 
with the above prediction. 
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FIG. 4: Specific heat versus temperature for a graph with a degree distribution that is Poisson (Q « 4). The plots refer to an 
average over 500 realisations for N — 10 3 , 10 4 and 50 realisations for N = 10 5 . 
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FIG. 5: Scatter plot of single node magnetisations versus degree for one network of iV 
temperatures: (a) below T c (k B T c /J = 2.87 ± 0.01), (b) above T c . 
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IV. THE CRITICAL TEMPERATURE 



In this section we discuss two approaches, the Bethe-Peierls approximation and the replica method, that allow us to 
get precise results for the critical temperature. Both approaches assume that there are no degree correlations present 
in the network. 



(a) (b) 



E 




FIG. 6: Scatter plot of magnetisations of all nodes of degree k = 2 for the same network as in Fig. 03 versus node index (i.e. 
the time step at which this node is added when the network grows according to the BA-rule). Diamonds represent those nodes 
connected to hubs whereas crosses represent nodes of degree 2 not connected to hubs. Older nodes (smaller index) are on 
average more connected to hubs than later ones. Different clusters represent different temperatures. In (a) we show the results 
for two temperatures below T c while in (b) the temperatures are above the critical one. The inversion phenomenon described 
in the text can clearly be seen. 
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FIG. 7: The specific heat (left) and the magnetisation (right) as a function of temperature for a BA network with /j, — 1/3. 
The average connectivity is Q = 4. 



A. Bethe-Peierls approximation 

The Bethe-Peierls (BP) approximation is the simplest of the cluster variation methods (2^|. It amounts to approx- 
imating the entropic part of the free energy, F(f3), by restricting the probability distribution p of a configuration of 
N spins to a combination of single-site and nearest-neighbour pair distributions. For a particular fixed network, this 
leads to the following expression for F(j3) 

k ■ 

F(p)f3 « /3(h) + - h) pi 1 * pi + \ EE E pa ln Pv ( 23 ) 



Sj=±l i j — 1 Si,8j=±l 



where 

H 



1 hi 

^EE-V"*' (24) 



i j=l 



Pi = 2 ( ' 

/'..- = ^ " " (26) 

Here = (siSj). Expressions for the local magnetisations mj and the neighbour spin-spin correlations ipy are 
obtained by looking for extrema of the free energy 

dF 

V* (27) 



dm, 
dF 



Vi,j (28) 



It can be shown that {27J| and |(2HJ| are fully equivalent to and J20J) . 

Next, the resulting selfconsistency equations are linearised in the m.;. In this way one obtains for the spin-spin 
correlations 

^ =tanh(Jij/fc B T) (29) 
After summing over all the vertices, the equation for the magnetisation becomes 

1 ]T(fc r - IK = ^E m -E TT^r- ( 3 °) 

Since we now want to get a simpler analytic expression for the critical temperature, we do not solve (|30H on single 
graphs, as already done in a more complete setting with the cavity approach. Instead, we average it over the network 
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realizations. This we do in two steps. Firstly, using a similar approximation as in section 2.1, we replace the second 
sum on the right hand side by a sum over all nodes. Hence, we obtain 

1 1 N ( 1 \ 

- - l)m r = - £ m r g> rj J ( 31 ) 

where again p r j = k r kj/(QN). Secondly, we insert the proportionality between rrij and fcj ^ implied by IT3j l: 

uij oc fcj ^- In the limit iV — > oo the resulting equation can again be written in terms of the average over the degree 
distribution. We then finally obtain 

oo OO OC . 

Q p (k)(k - l)* 1 "" = P(ki)kl» P ^ f JQ2 , \ ( 32 ) 

k=m fc 1= m fc 2 =m 1 + tanh [ kBTk » k ») 

For the case /x = 0, and for Q large, can be approximated and gives the estimate fcsT c /J w (Q2 — Q)/Q, m 
agreement with exact results [13] • Also, for the SAN on a BA-network, we get for Q large: fcsT c /J On a 

regular lattice with coordination number Q, the Bcthe-Peierls approximation gives J/(ksT c ) = In J /2 H3- As 

can be seen in table 1, this approximation gives a very good value for T c for Q > 4. 

Equation l|32(l can be solved numerically. In table 1 we present our results for the SAN on a BA network with 
P(k) = 2m(m + l)/k{k + l)(k + 2) (for this case, m = Q/2). 



TABLE I: Critical temperature fcsT c /J as a function of Q for the SAN on networks, as obtained from various approaches. 
On a regular lattice with coordination number Q, the Bethe-Peierls (BP) approximation gives J/fcsT c = In (Q/(Q — 2)) /2. 
This value is given in the second column. In the third column the results of the cavity approach are given, while the fourth 
column presents the estimate obtained within the Bethe-Peierls approximation on the network, equation do 2 1 . The fifth column 
gives the solution to 1561 obtained in the replica approach. The last column gives the estimates coming from the Monte 
Carlo simulations. The results for the cavity method and the simulations were obtained on BA networks. The results for the 
network BP and the replica method were obtained on uncorrelated networks. This has dramatic consequences for Q — 2. For 
Q — 2, T c = on BA networks because BA networks with Q — 2 are simple trees without loops. In contrast, uncorrelated 
networks with Q — 2 may feature loops as well as disconnected parts. 



Q lattice BP 


cavity 


network BP 


replica 


Monte Carlo 


2 


0. 


0. 


0.95230 


0.94614 


0. 


4 


2.8854 


2.87T 0.01 


2.95508 


2.92468 


2.91 ± 0.02 


6 


4.9326 


4.92± 0.01 


4.9554 


4.94511 
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6.9521 


6.94± 0.01 


6.9406 


6.95796 




10 


8.9628 


8.94± 0.01 


8.9087 


8.96607 


8.96 ± 0.01 


20 


18.9824 


18.98 ± 0.01 


18.9107 


18.98187 


19.09 ± 0.04 



B. Replica theory for an uncorrelated network 



The replica approach is a powerful mathematical technique that provides a way to perform an average over random 
network ensembles. For the present model, and in absence of degree correlations, it leads to a relatively simple 
analytical equation from which the critical temperature can be determined numerically. 

In order to define the random ensemble of networks we start by defining the adjacency matrix C of a graph. This 
matrix is the N x N matrix whose element cy is one if there is a link between node i and j and zero otherwise. 
Clearly, 

JV 

^ ('ij = ki (33) 
3=1 

Moreover, for undirected graphs C is symmetric, cy = Cjj. Using the adjacency matrix we can write the Hamiltonian 
TL for our model as 

j N N 

H = ~- ^2 -I'l^k,.!,- , )<;,■%■■< , - H^Si (34) 

i^tj i=l 
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where in our case $(fej, kj) = Jij/J = Q 2/i / '(kikj)^ . The replica approach can, however, be performed for more general 
forms of In Eq. (|34[) . we have also included an external magnetic field H for later convenience. 

We calculate the typical properties of all networks with a given fixed set of degrees {ki}. The matrix elements c^- 
are assumed to be completely uncorrelated apart from the constraint i|33|) . Therefore given that 



1 I „ 



we can write their joint probability distribution as 

men}) oc n Pi**) n 5 1 e <=« - fc * 



(35) 



(36) 



As usual in disordered systems, we compute the quenched average of the free energy density per site, f q (/3), from 
which typical properties can be determined 



1 



(3f q (f3) = lim -[InZ] 



(37) 



where (3 = l/(fcsT). Since the calculation is rather involved but standard [bl l2ll l30|. we summarize here only the 
major steps and results. Using the equality InZ = lim„^o h(Z n ~ 1) the quenched free energy density is written as 



f3f g ((3)= lim hm-L([Z" (/?)]-!) 



The replicated partition function Z n ((3) equals 

r /3 



n JV 



EE J ^ Q +/^EE s " 



a— 1 i,j 



a— 1 i= 1 



(38) 



(39) 



where Sj = {s*, . . . , s™ }. For the distribution of connectivities given in Eq. Q36[) . the average [•] is 





N 




Y[dcijP(dj) 




E c ' j ^ 




i=l 





where A/" is the normalisation 



N = 





N 




\{dc l3 P{a l3 ) 


IP 


E c ' j ^ 


i<j 


i=l 





In order to calculate this average, it is common to introduce an exponential representation of the constraint 



N 

S ( °ij 



2tt 



(40) 



(41) 



(42) 



A straightforward but lengthy calculation in which we integrate over the disorder and the auxiliary variables allows 
us to obtain the free energy density in terms of the functional order parameter 



1 N 



(43) 



and its canonical conjugate R k - This order parameter is an order parameter in the replica space which is the joint 
density of finding a spin configuration s with average connectivity k at each site, when a link has been removed from 
that site. The result is 



Pfq(P) - Extr {R k (s),R k (s)} 



£ P(k) m]T [R k m V*^ + f£E R k (s)R k , (a 



a /3J®(k,k')s-a 



(44) 



k,k' s,<7 
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where we also performed the rescaling iR k (s) — > —QRk{s)- 

Stationarity of the free energy with respect to Rk(s) and Rk(s) leads to the saddle-point equations 



R k (s) 



kP(k) 


R k (s) 


k- 


1 e PHT,l = ^ a 


Q E, 


Rk(<?) 





(45) 



If all the <&(fc, k') are positive, so that the model only has ferromagnetic interactions, it can be expected that the 
solution to these equations has replica symmetry (RS). For the order parameter and its conjugate, this assumption 
takes the form 



R k (s)= / dhW k (h) 



and 



Rk{s) = / duQ k (u) 



[2cosh03/i)f 



[2 cosh(/3u)] n 



(46) 



(47) 



respectively. Here /i and u are the local cavity message and the propagated field that we already encountered in 
section 2. While the cavity approach calculates these fields for a specific network, in the replica approach we can 
obtain their probability distributions over the set of all networks obeying the constraints (|33|) . These distributions 
are denoted as W k {h) and Q k {u) respectively. Within the RS ansatz, the saddle-point equations become 



Qk(u) = J dhW k '(h)S \u- itanlT 1 [tanh(/3/i) tanh [/?J$(fc, fc')]] 



W k {h) 



P{k)k 



Q 



'fc-i 



fc-i 



S ih-J^ui- H 



1=1 



duiQ k (ui) 

and the free energy density equals 

Pfqlfi) = -Qln2 + Q^ I dudhQ k {u)W k {h)\n[l + tanh(/3u) tanh(/3/i)] 



(48) 
(49) 





" k 




/ 


Y\_duiQ k (ui) 


In 









fe 

/ dhdtiW k {h)W k >{ti)Fkk>{P,h,ti) 



2 cosh 



" =1 ui + H 



Ili = i 2 cosh (/?«,) 



(50) 



where 

F kk > (P, h, ti) = In [cosh [/3J$(fc, fc')] + sinh [f3 J$(fc, fe')] tanh(/?/i) tanh(/%')] 
The average zero-field magnetisation per site, M, then follows immediately 



M 



dH 



(Jf = 0) = ^P(fc) | 



J^duiQ fe ( 



i=i 



tanh ( P ^ 



Ui 



(51) 



1=1 



From Eq. (|48|) and Eq. I|49|l we can obtain an equation that contains only the functions Qk(i 



^ P(q)q f 



9-1 

, [ duiQ q (ui) 
.1=1 



6(u~G kq (p)) 



(52) 



12 



where 

i . r / «zi \ i 

(53) 



G kg (P) = -tann" 1 



tanh ( P^ui J tanh [/3 J$(fc, q)\ 



i=i 



One can try to solve this set of equations, for example using population dynamics techniques. Here we only 
investigate the simpler question of locating the critical temperature. For this, we assume that the cavity fields u are 
very small near the transition. The argument of the delta-function in Eq. (|52[1 can then be linearised using 

9-1 

G kq (J3) «tanh [/?J$(fc, q)\ ^ m 

i=i 

We next multiply both sides in Eq. (|52|) by u and we integrate. This gives 

9-1 



/ duQk(u)u = J2^^ taa MPJ®(k,q)]Yl I dui ® 

J q Q 1=1 J 



(ui)u t (54) 



Finally, we denote Xk — J Qk{u)udu and define a matrix A(f3) with elements 



A q m(P) = P ^ q ^ — ^ tanh [/3J$(fc, q)] q,k>m (55) 

(m is the smallest degree appearing in the network). Finding the critical temperature now amounts to locating the 
value of (3 for which the matrix A has an eigenvalue 1, i.e. to solving 

det (A{p c ) -I) = (56) 

In principle, the matrix A is infinite dimensional. For a given P(k),m and $(fc,/c') one can, however, calculate an 
estimate T C (K) for the critical temperature by truncating the matrix and limiting q and k to be smaller than a given 
K. An extrapolation for K — > oo then gives T c . We have performed such a calculation for the SAN for different values 
of Q. Some numerical values can be found in table 1. 



V. MONTE CARLO SIMULATION RESULTS 



Finally, we also investigated our model numerically on a BA network. This allows us to investigate effects that are 
neglected within the cavity approach, such as the appearance of loops in the network. 

Firstly, we investigated the SAN with various values of Q. When a BA network is grown by the addition of m new 
links, one always ends up with an even value of Q, since Q = 2m. In order to obtain an odd average connectivity, 
alternating values of m were used. For example, by alternating m between 1 and 2 we obtain Q = 3. We used the 
standard Metropolis algorithm to simulate our model. Typically, we averaged over 400 uncorrelated spin configurations 
at each temperature and for each realisation of the network. 

The simulations were made for networks with N ranging between 100 and 56000. At low values of Q we found large 
fluctuations in the properties of different realisations. We typically investigated between 200 and 5000 realisations of 
the network, depending on the size of the network and the quantity under investigation. 

The value of T c was obtained from the intersection of the Binder cumulants Un = 1 — (Af 4 )/(3(M 2 ) 2 ), with M 
the total magnetisation in a given configuration. In order to obtain a value for the intersection point, we first fitted 
a curve through the data points (generally a fifth degree polynomial) and then solved for the intersections of the 
fitted curves. From this, we obtain the finite-size critical temperature' T C (N). These were then extrapolated using 
the standard relation 

T c (A0 ~ T c (oo) - 6iV-2^( 1 +^) (57) 

where w is a correction-to-scaling exponent, introduced by Binder [3ll l32l l33l | . and v the "correlation length" exponent 
proper to lattice models. In our network models the product wv figures as a single number, since a correlation length 
is not defined. In Fig. |SJ a typical example of such a fit is shown. The resulting estimate of T c , together with those 
obtained for other Q-values, is given in table 1. Moreover, from the data in Fig. |H1 the estimate (1 + ojv)/(2 — a) = 
0.58 ± 0.12 is obtained. The network with Q = 3 was only investigated with the simulation method and for this case 
we find k B T c /J = 1.85 ± 0.03. 
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FIG. 8: Finite-size estimate of the critical temperature (divided by the mean field estimate JQ/fcs) as a function of 1/N 
obtained using the intersections of the cumulants. An example of such an intersection is given in the second plot. The average 
connectivity is Q — 4. 



From the Binder cumulant, we can estimate the exponent a. Indeed, the derivative of the cumulant U' N — dU^ /dT 
at T C (N) scales as 



U'n(Tc) 
U' N ,(T C ) 



(58) 



The derivative can be obtained from the polynomial fit used to determine T C (N). The values that we obtain are 
consistent with the assumption that a — 0. 

Next, we calculated the magnetisation at T C (N) as a function of N. This allows us to estimate the exponent 
P/(2 — a) as ~ 0.25. In Fig. [5] (left side) we show our numerical results for Q = 10. This result is consistent with 
the mean field values (3 — 1/2 and a = 0. As a further test, we plotted the squared magnetisation as a function of 
T . Below T c , we expect this to be a linear function and from a fit to this form we can get a second estimate of T c 
(see Fig. [51 right side). The values that we find in this way are in agreement with those coming from the cumulants, 
although the precision is less good. 




4.0, 4.5 

Log 10 (N) 




0.4 0.6 0.8 

T/T(MF) 



FIG. 9: iV-dependence of the magnetisation of the SAN at T C (N) (Q = 10). The slope is -0.25 ± 0.01. The right side shows 
the magnetization squared as a function of temperature for a network of 3200 nodes, averaged over 400 samples. 

We also computed the finite-lattice susceptibility \' i introduced by Binder (jyj 

X ' cx (Af 2 ) - (|M|) 2 (59) 

By taking the modulus of the magnetisation, one tries to minimise finite-size effects. Besides giving information 
about T c , which is in agreement with the results coming from the cumulants and the magnetisation, the susceptibility 
provides information on the critical exponent j s . At T c , one expects \' cx jV 7s /( 2-Q ). We calculated x' as a function 
of N at a temperature close to T c . We find that 7 s /(2 — a) is close to 0.47 ± 0.03, which is consistent with j s = 1 
using the earlier estimate of a. 

The plots for the specific heat as a function of temperature (Fig. left side) look similar to those obtained from 
the cavity method. The specific heat does not diverge with size but saturates, which implies a ~ 0. However, the 
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plots do not show the usual, mean field jump at T c , but are in agreement with the presence of logarithmic corrections. 
The maximum in the specific heat below T c that was found within the cavity approach (Fig. [3J is also clearly visible 
in the data. When we change the value of /i from 1/2 to 1 we do find results consistent with the usual mean field 
behaviour, i.e. with a jump at T c (see Fig. 10, right side and compare with Fig. 4). Indeed fi = 1 corresponds to 
7' — > 00, which corresponds to a network with a narrow degree distribution. 



1.6- 
1.2- 
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T/T (MF) 
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o 
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FIG. 10: The specific heat as a function of temperature for a network of 5600 nodes (SAN model , left) and for the modified 
model with y, = 1 (right). Averages were taken over 400 samples. The average connectivity is Q — 10. 

We also performed a completely similar investigation of the network with fi = 1/3 and Q = 10. In Fig. 1111 we show 
logdog plots of the magnetisation (left) and the susceptibility (right) as a function of N at the numerically determined 
value of T c . From these data we obtain /3/(2 - a) = 0.358 ± 0.001 and j a /(2 - a) = 0.282 ± 0.002. These results arc 
consistent with the predictions coming from 112(1 : (3/ (2 — a) — 1/3 and 7 s /(2 — a) = 1/3. The remaining discrepancy 
is probably due to errors in the determination of T c . We also obtained data for the specific heat, but they do not allow 
a precise determination of a. From 1(12(1 . the prediction a = — 1 follows. This implies that the specific heat decreases 
linearly to its high-temperature background value. The data shown in Fig. 1121 are consistent with this expectation. 
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FIG. 11: Log-log plot of the magnetisation (left) and the susceptibility (right) versus N at T c for the case fx = 1/3, Q — 10. 

From the simple mean field theory of section 2, it follows that the local order parameter m, oc We also used 

this relation in our derivation of the BP approximation. In order to get further confirmation for this type of scaling, 
we checked it numerically. Our results are shown in Fig. 1131 for the cases fi = 0, 0.5 and 1 on a BA-network. For 
high values of the degree, the statistics is not very good, but for smaller fc^-values the scaling 1(14(1 seems to be well 
satisfied. 

We finally remark that for the SAN with Q = 2, the Monte Carlo simulations indicate a zero critical temperature. 
Indeed, for this value of Q, the BA network has the structure of a tree and therefore it cannot support order. This 
is in contrast with the results coming from the replica approach. However, in that approach loops are not excluded 
and the network can be a collection of disconnected clusters, the largest of which determines the non-zero T c in the 
large- TV limit. To test this, we generated a set of random (uncorrelated) networks with Q = 2 and 7 = 3 by ascribing 
to each vertex a degree taken from the distribution P(k). The links going out of the vertices are then randomly paired 
to create a network. The resulting network is indeed found to consist of disconnected parts. In Fig. ^]we show 
Monte Carlo results for the Binder cumulant as a function of temperature for two networks (N — 2000 and TV = 4000) 
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FIG. 12: The specific heat as a function of temperature for a network of 1000 nodes (fi — 1/3). The average connectivity is 
Q = 10. Notice the linear regime for 8 < k B T / J < 10.2 




FIG. 13: The magnetization of individual spins as a function of degree for a network of 1000 nodes and different values of /i. 
The averages are made over 1000 networks. The values of fi and the predicted fc-dependence, respectively, is: /j, = - linear; 
fi — 0.5 - square root; p = 1 - constant. The temperature is slightly below the respective critical temperatures. 



generated in this way. We took /x = 1/2. Clearly there is an intersection from which we obtain the finite-size estimate 
for the critical temperature. After extrapolation (TV — > oo) we obtain from these kind of data fcsT c /J = 0.92648. For 
comparison, the BP prediction for kT c /J is 0.952 and the replica method gives 0.946 (Table 1). 




0.450 0.455 0.460 0.465 0470 0.475 

T/T (MF) 



FIG. 14: Binder cumulants as a function of temperature for a network with 7 = 3, Q = 2 and fj, = 1/2. The figure shows a 
clear intersection of the cumulants for iV = 2000 (circles) and N = 4000 (squares) indicating the presence of a critical point at 
finite temperature. The averages were done over 6000 and 5000 realizations of the networks, respectively. The fit with a third 
degree polynomial is also shown in the figure. 
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VI. NONEQUILIBRIUM MODELS 

In this section, we show that the basic result of this paper, (|12f> . is also true for a noncquilibrium model. In particular, 
we will investigate the contact process [35j | on a network. This process is a well known model originally introduced to 
describe the spreading of a disease in a population. Later, it was found to be related to directed percolation (3(| and 
became the prototype model used in the study of absorbing state phase transitions [13 . Recently the contact process 
has also been applied in metapopulation ecology [381] . but in this paper we will use the epidemiological language. 

On a network we define the contact process as follows. Each vertex i can be in two states that we denote as ill 
(rii = 1) and healthy (n,; = 0). The dynamics of the model is given by a continuous time Markov process [3fj. An ill 
site can cure with a rate 1 (this fixes the time scale). A healthy site becomes ill with a rate that equals A times the 
number of ill neighbours, where on a network two vertices are called neighbours' if they are linked. 

This model has an absorbing state: if all sites are healthy, they will stay healthy forever. For an infinite system, 
there is a phase transition at some A c . If A < A c , the system will evolve to the absorbing state. For A > A c , there 
will always be a finite density of ill' sites. This density is the order parameter of the model and it goes to zero if 
one approaches A c from above. Vespignani and coworkers showed that on a scale-free network for 2 < 7 < 3, 

A c = and the exponents are 7-dependent. For 3 < 7 < 4, A c > and the critical exponents are again 7-dependent. 
Finally, for 7 > 4, A c > and critical exponents assume the mean-field values for the contact process on a regular 
lattice. This scenario is reminiscent of that for the opinion formation model. 

We next generalise the contact process by taking the infection rate as 

Ay = \Q 2 »(hk 3 )-» (60) 

Using standard techniques from the theory of stochastic processes [39|. one can show that pi = {rii) (where (.) in this 
paragraph denotes the average over histories of the stochastic process) obeys the exact equation 

$--^+E A ^( 1 -"^) ( 61 ) 

The sum runs over the neighbours of i. We now make two approximations similar to those performed in section 2. 
The first one, the dynamical mean field approximation, puts ((1 — rii)ni) w ((1 — ni))(ni) = (1 — Pijpi- Next, again 
following Bianconi [lj|, we replace the sum over the neighbours by a sum over all the nodes, and take for Ay its 
average over all realisations of the network with a given set of degrees. Then (|61f) becomes 

d N 

-j± = -p l + {l-p l )Y J [\ l] ]p J (62) 
3=1 

Now (see section 2) [Ay] = Ay-py with py = kikj/(QN), the probability that nodes i and j are connected. This gives 

Pi + d-pd^^kl-^k]-^ (63) 



dt rt v N 

We introduce O 



e =V5>r^ (64) 

From here on, we will assume that pi depends only on the degree of i. We can then write for a large network 

Q = Q^- 1 Y,P{k)k 1 -»h (65) 

where pk is the density of sites that are ill and have degree k. 
Using (|64|l and this assumption, l|63|l becomes 

^ = -p fe + A(l-p fe )Q"fc 1 ^e (66) 
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We are interested in the static properties, i.e. the properties of the model in the long time limit. Then, from 166(1 
we get 

Pk - i + AQMfci-,e (67) 

Using l|65|l . we find a self-consistency equation for 

fe 



e = Ag 2 ^- i yp(fc)- — L- -f - (68) 



We can now in principle analyse this equation. 

But as was the case for the opinion formation model, it is easier to transform away the /i-dependence. This is most 
easily done by rewriting (|68|l as an integral 

r°° 1.2-2^0 

@ = AXQ » L r1 ITA^e* ^ 

where A is a normalisation constant and where the power-law form for P(k) is inserted. If we change variables to 
k' = Q^k 1 -^, JTHJ becomes (if fi < 1) 

e = A 'L (k '^TTWe dk ' <™> 

where A' is another constant and m! = m 1_M Q M . 

This is precisely the form l(69|) assumes for // = provided we shift 7 to 7' given by 

7 = I-77 (71) 

Thus, the same relation holds as for the static Ising case But notice that now, taking fi = 1/2 and a Barabasi- 
Albert network (i.e. the epidemiological version of the SAN), leading to 7' = 5, we are in a situation where mean-field 
theory should hold without logarithmic corrections. 

The relation (|12f> thus appears to be quite general and one may wonder whether there are any exceptions to it. We 
first discuss a recently proposed realization of the Bak-Tang-Wiesenfeld sandpile model on scale-free networks 43] . 
In this model an avalanche can be generated through the following dynamical rules: 1) at each time step, a grain is 
added at a randomly chosen network node i: 2) if the height hi at node i exceeds a given threshold Zi, it becomes 
unstable and an integer number of grains, n[zi], at the node topple, with n[zi] — 1 < Zi < n[zi], to randomly chosen 
n[zi] nodes among fcj adjacent ones; 3) whenever adjacent nodes become unstable toppling takes place also there, on 
all unstable nodes in parallel, and the avalanche continues until there are no unstable nodes left. In the proposed 
version of the model [43| , a parameter 77 is introduced and the threshold of node i is taken to be 

Zi = (72) 

Note that for r\ = the model features a simple degree-dependent threshold, and in the extended model k is replaced 
by fc 1_?) . Note that this is reminiscent of the transformation k — > which we discussed in other contexts in this 
paper. 

A key quantity in the description of the avalanche dynamics by mapping each avalanche to a tree, is the branching 
probability P(,(n]_ that a node, which receives a grain from a neighbour, generates n branches. This probability consists 
of two factors [iij. If branching occurs for a given node i, n[zi] branches are generated. Therefore, the first factor, 
Pi(n), is the probability that n coincides with n[zi] for a node i already connected to the tree. The second factor is 
the probability that the height hi takes the value n — 1 at the moment that node i receives a grain from a neighbour. 
This probability, P2(n), is not important for us here, but in the case of independent random heights 0, 1, n — 1 in 
the inactive state of the sandpile, it equals 1/n. If we would adopt a continuum approximation, treating n as a real 
number, and approximate n[zi] by zi, we would arrive at the result 

P h (n) = (n 1 /^P(n 1 /^)/Q)P 2 (n), (73) 

where the first factor, Pi(n) = kP(k)/Q, is the probability that the first neighbour of a node has degree fc, with in 
our case k = n 1 ^ 1 ~ ,1 \ This probability differs from P(k) in that it presupposes that the generating node is already 
connected to the tree. 
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If we compare this Pb{n) with its counterpart for the model with r\ = 0, we easily observe that the extended model 
is equivalent to the basic model with rjf = 0, provided the topological exponent is transformed according to 

7 ' = T~^7j < 74 > 

This is in agreement with our general exponent relation, but, as we shall show now, and as Goh et al. already found 
p3| this result is incorrect. 

The correct exponent relation is found when taking into account properly that n is always an integer, and that all 
thresholds in the interval n — 1 < Zj < n contribute to the probability of generating n branches. P\(n) thus consists 
of a sum, or, for simplicity and without loss of relevant precision, an integral, so that 

Pb{n) = ( / kP(k)dk/Q)P 2 {n), (75) 

J(n-l)VCi-ti) 

from which follows the correct relation 

, 7-2?? , . 

7 = 1-7' (76) 

which was already obtained by Goh et al |43| . We conclude that an exception to our general exponent relation arises 
here due to the discrete (integer) character of the toppling process. Indeed, the general relation is recovered in a 
(rough) continuum approximation of the problem. 

A second exception to our general exponent relation is provided by the study of Dezso and Barabasi 01 01 disease 
spreading on a scale-free network. The model they consider starts from the usual contact process, for which it is known 
that the epidemic threshold A c vanishes for 7 < 3. This is similar to the divergence of the critical temperature for an 
Ising model with constant interactions on a scale-free network. We have seen that topology-dependent interactions 
are a way to get around this and we have discussed this in detail also for the contact process. 

However, Dezso and Barabasi provide an alternative means of rendering the epidemic threshold finite and thus 
offering new avenues for controlling diseases (e.g., eradicating viruses). Their strategy consists of curing the hubs with 
a probability that scales with the degree A; of a node as k a . Note that a = corresponds to the usual model in which 
all nodes are cured with the same probability. 

As a result Eq. (|69|) takes the modified form 0] 

r°° k 2 f) 

e = AXQ L fc_7 F+A*e*' (77) 

from which we can derive the following exponent relation 

y = 7 -^, (78) 

1 — a 

with 7' the effective topological exponent for an equivalent model with degree-independent curing rate (a' = 0). 
Curiously, this new relation is akin to that of the previous exception (sandpile model), but this coincidence has to 
our insight no significance. 



VII. DISCUSSION 



In this paper, we investigated the critical properties of an Ising model and a contact process with topology-dependent 
interactions on scale-free networks. The interaction strength between two spins, or the infection rate between two 
individuals, was assumed to be proportional to (fejfej) - '*. We have developed mean-field theories for these models from 
which our main result follows: the critical behaviour can always be related to that of the corresponding model with 
homogeneous couplings (p/ = 0) on a network with a modified degree distribution 7', where 7' is given by equation 
(|12fl . Due to the small-world property of scale-free networks, mean field theory is generally believed to be exact. This 
expectation is found to be true also in the present case. We performed extensive numerical calculations for the Ising 
model on a BA-network, focusing on the cases fi = 1/2 (the SAN, which was the original motivation for the present 
work) and /1 = 1/3. We used two techniques: numerically exact' studies with the cavity approach, and Monte Carlo 
simulations. The first approach assumes the absence of loops in the network, which is a good approximation for a 
BA-network. The Monte Carlo calculations approximate the true thermal average by an average over a finite set 
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of well chosen spin configurations. Within their numerical accuracy, both approaches give the same results. More 
importantly, they provide strong evidence that the exponent relation 112fl is correct. 

For static network ensembles it is possible to obtain an exact equation for the location of the critical temperature 
from the replica approach, see equations (|55|l and H5(jf) . A simpler equation for this quantity, l|32[l . can be derived from 
the Bethe-Peierls approximation. From the numerical values listed in table 1, it can be seen that this simple aproach 
gives results that are accurate to better than one percent. From the table one also observes that the differences in 
the critical temperature as obtained from the different approaches are very small for Q > 4. 

Besides the contact process, the behaviour of several other (non) equilibrium models have been studied on scale-free 
networks. In particular, we mention here diffusion-annihilation |45j | and the Bak-Sneppen model |4rj| . In both cases, 
the critical behaviour shows a 7-dependence that is reminiscent of that of the Ising model and the contact process. It 
remains to be investigated whether appropriately extended versions of these models obey the relation (|12l) or whether 
they follow modified transformation laws as we found to be the case for the sandpilc model of ref. 43 or the modified 
contact process of ref. 44. 
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